home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 11 / Info-Mac_XI_Disc_1.cdr_ / Info-Mac XI Disc 1.cdr / Programs / Science & Math / MacEspresso 1.0 / espresso / matrix.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-12-10  |  11.4 KB  |  510 lines  |  [TEXT/R*ch]
  1. #include "port.h"
  2. #include "sparse_int.h"
  3.  
  4. /*
  5.  *  free-lists are only used if 'FAST_AND_LOOSE' is set; this is because
  6.  *  we lose the debugging capability of libmm_t which trashes objects when
  7.  *  they are free'd.  However, FAST_AND_LOOSE is much faster if matrices
  8.  *  are created and freed frequently.
  9.  */
  10.  
  11. #ifdef FAST_AND_LOOSE
  12. sm_element *sm_element_freelist;
  13. sm_row *sm_row_freelist;
  14. sm_col *sm_col_freelist;
  15. #endif
  16.  
  17.  
  18. sm_matrix *sm_alloc(void)
  19. {
  20.     register sm_matrix *A;
  21.  
  22.     A = ALLOC(sm_matrix, 1);
  23.     A->rows = NIL(sm_row *);
  24.     A->cols = NIL(sm_col *);
  25.     A->nrows = A->ncols = 0;
  26.     A->rows_size = A->cols_size = 0;
  27.     A->first_row = A->last_row = NIL(sm_row);
  28.     A->first_col = A->last_col = NIL(sm_col);
  29.     A->user_word = NIL(char);        /* for our user ... */
  30.     return A;
  31. }
  32.  
  33.  
  34. sm_matrix *sm_alloc_size(int row, int col)
  35. {
  36.     register sm_matrix *A;
  37.  
  38.     A = sm_alloc();
  39.     sm_resize(A, row, col);
  40.     return A;
  41. }
  42.  
  43.  
  44. void sm_free(sm_matrix *A)
  45. {
  46. #ifdef FAST_AND_LOOSE
  47.     register sm_row *prow;
  48.  
  49.     if (A->first_row != 0) {
  50.     for(prow = A->first_row; prow != 0; prow = prow->next_row) {
  51.         /* add the elements to the free list of elements */
  52.         prow->last_col->next_col = sm_element_freelist;
  53.         sm_element_freelist = prow->first_col;
  54.     }
  55.  
  56.     /* Add the linked list of rows to the row-free-list */
  57.     A->last_row->next_row = sm_row_freelist;
  58.     sm_row_freelist = A->first_row;
  59.  
  60.     /* Add the linked list of cols to the col-free-list */
  61.     A->last_col->next_col = sm_col_freelist;
  62.     sm_col_freelist = A->first_col;
  63.     }
  64. #else
  65.     register sm_row *prow, *pnext_row;
  66.     register sm_col *pcol, *pnext_col;
  67.  
  68.     for(prow = A->first_row; prow != 0; prow = pnext_row) {
  69.     pnext_row = prow->next_row;
  70.     sm_row_free(prow);
  71.     }
  72.     for(pcol = A->first_col; pcol != 0; pcol = pnext_col) {
  73.     pnext_col = pcol->next_col;
  74.     pcol->first_row = pcol->last_row = NIL(sm_element);
  75.     sm_col_free(pcol);
  76.     }
  77. #endif
  78.  
  79.     /* Free the arrays to map row/col numbers into pointers */
  80.     FREE(A->rows);
  81.     FREE(A->cols);
  82.     FREE(A);
  83. }
  84.  
  85.  
  86. sm_matrix *sm_dup(sm_matrix *A)
  87. {
  88.     register sm_row *prow;
  89.     register sm_element *p;
  90.     register sm_matrix *B;
  91.  
  92.     B = sm_alloc();
  93.     if (A->last_row != 0) {
  94.     sm_resize(B, A->last_row->row_num, A->last_col->col_num);
  95.     for(prow = A->first_row; prow != 0; prow = prow->next_row) {
  96.         for(p = prow->first_col; p != 0; p = p->next_col) {
  97.         (void) sm_insert(B, p->row_num, p->col_num);
  98.         }
  99.     }
  100.     }
  101.     return B;
  102. }
  103.  
  104.  
  105. void sm_resize(register sm_matrix *A, int row, int col)
  106. {
  107.     register int i, new_size;
  108.  
  109.     if (row >= A->rows_size) {
  110.     new_size = MAX(A->rows_size*2, row+1);
  111.     A->rows = REALLOC(sm_row *, A->rows, new_size);
  112.     for(i = A->rows_size; i < new_size; i++) {
  113.         A->rows[i] = NIL(sm_row);
  114.     }
  115.     A->rows_size = new_size;
  116.     }
  117.  
  118.     if (col >= A->cols_size) {
  119.     new_size = MAX(A->cols_size*2, col+1);
  120.     A->cols = REALLOC(sm_col *, A->cols, new_size);
  121.     for(i = A->cols_size; i < new_size; i++) {
  122.         A->cols[i] = NIL(sm_col);
  123.     }
  124.     A->cols_size = new_size;
  125.     }
  126. }
  127.  
  128.  
  129. /*  
  130.  *  insert -- insert a value into the matrix
  131.  */
  132. sm_element *sm_insert(register sm_matrix *A, register int row, register int col)
  133. {
  134.     register sm_row *prow;
  135.     register sm_col *pcol;
  136.     register sm_element *element;
  137.     sm_element *save_element;
  138.  
  139.     if (row >= A->rows_size || col >= A->cols_size) {
  140.     sm_resize(A, row, col);
  141.     }
  142.  
  143.     prow = A->rows[row];
  144.     if (prow == NIL(sm_row)) {
  145.     prow = A->rows[row] = sm_row_alloc();
  146.     prow->row_num = row;
  147.     sorted_insert(sm_row, A->first_row, A->last_row, A->nrows, 
  148.             next_row, prev_row, row_num, row, prow);
  149.     }
  150.  
  151.     pcol = A->cols[col];
  152.     if (pcol == NIL(sm_col)) {
  153.     pcol = A->cols[col] = sm_col_alloc();
  154.     pcol->col_num = col;
  155.     sorted_insert(sm_col, A->first_col, A->last_col, A->ncols, 
  156.             next_col, prev_col, col_num, col, pcol);
  157.     }
  158.  
  159.     /* get a new item, save its address */
  160.     sm_element_alloc(element);
  161.     save_element = element;
  162.  
  163.     /* insert it into the row list */
  164.     sorted_insert(sm_element, prow->first_col, prow->last_col, 
  165.         prow->length, next_col, prev_col, col_num, col, element);
  166.  
  167.     /* if it was used, also insert it into the column list */
  168.     if (element == save_element) {
  169.     sorted_insert(sm_element, pcol->first_row, pcol->last_row, 
  170.         pcol->length, next_row, prev_row, row_num, row, element);
  171.     } else {
  172.     /* otherwise, it was already in matrix -- free element we allocated */
  173.     sm_element_free(save_element);
  174.     }
  175.     return element;
  176. }
  177.  
  178.  
  179. sm_element *sm_find(sm_matrix *A, int rownum, int colnum)
  180. {
  181.     sm_row *prow;
  182.     sm_col *pcol;
  183.  
  184.     prow = sm_get_row(A, rownum);
  185.     if (prow == NIL(sm_row)) {
  186.     return NIL(sm_element);
  187.     } else {
  188.     pcol = sm_get_col(A, colnum);
  189.     if (pcol == NIL(sm_col)) {
  190.         return NIL(sm_element);
  191.     }
  192.     if (prow->length < pcol->length) {
  193.         return sm_row_find(prow, colnum);
  194.     } else {
  195.         return sm_col_find(pcol, rownum);
  196.     }
  197.     }
  198. }
  199.  
  200.  
  201. void sm_remove(sm_matrix *A, int rownum, int colnum)
  202. {
  203.     sm_remove_element(A, sm_find(A, rownum, colnum));
  204. }
  205.  
  206.  
  207.  
  208. void sm_remove_element(register sm_matrix *A, register sm_element *p)
  209. {
  210.     register sm_row *prow;
  211.     register sm_col *pcol;
  212.  
  213.     if (p == 0) return;
  214.  
  215.     /* Unlink the element from its row */
  216.     prow = sm_get_row(A, p->row_num);
  217.     dll_unlink(p, prow->first_col, prow->last_col, 
  218.             next_col, prev_col, prow->length);
  219.  
  220.     /* if no more elements in the row, discard the row header */
  221.     if (prow->first_col == NIL(sm_element)) {
  222.     sm_delrow(A, p->row_num);
  223.     }
  224.  
  225.     /* Unlink the element from its column */
  226.     pcol = sm_get_col(A, p->col_num);
  227.     dll_unlink(p, pcol->first_row, pcol->last_row, 
  228.             next_row, prev_row, pcol->length);
  229.  
  230.     /* if no more elements in the column, discard the column header */
  231.     if (pcol->first_row == NIL(sm_element)) {
  232.     sm_delcol(A, p->col_num);
  233.     }
  234.  
  235.     sm_element_free(p);
  236. }
  237.  
  238. void sm_delrow(sm_matrix *A, int i)
  239. {
  240.     register sm_element *p, *pnext;
  241.     sm_col *pcol;
  242.     sm_row *prow;
  243.  
  244.     prow = sm_get_row(A, i);
  245.     if (prow != NIL(sm_row)) {
  246.     /* walk across the row */
  247.     for(p = prow->first_col; p != 0; p = pnext) {
  248.         pnext = p->next_col;
  249.  
  250.         /* unlink the item from the column (and delete it) */
  251.         pcol = sm_get_col(A, p->col_num);
  252.         sm_col_remove_element(pcol, p);
  253.  
  254.         /* discard the column if it is now empty */
  255.         if (pcol->first_row == NIL(sm_element)) {
  256.         sm_delcol(A, pcol->col_num);
  257.         }
  258.     }
  259.  
  260.     /* discard the row -- we already threw away the elements */ 
  261.     A->rows[i] = NIL(sm_row);
  262.     dll_unlink(prow, A->first_row, A->last_row, 
  263.                 next_row, prev_row, A->nrows);
  264.     prow->first_col = prow->last_col = NIL(sm_element);
  265.     sm_row_free(prow);
  266.     }
  267. }
  268.  
  269. void sm_delcol(sm_matrix *A, int i)
  270. {
  271.     register sm_element *p, *pnext;
  272.     sm_row *prow;
  273.     sm_col *pcol;
  274.  
  275.     pcol = sm_get_col(A, i);
  276.     if (pcol != NIL(sm_col)) {
  277.     /* walk down the column */
  278.     for(p = pcol->first_row; p != 0; p = pnext) {
  279.         pnext = p->next_row;
  280.  
  281.         /* unlink the element from the row (and delete it) */
  282.         prow = sm_get_row(A, p->row_num);
  283.         sm_row_remove_element(prow, p);
  284.  
  285.         /* discard the row if it is now empty */
  286.         if (prow->first_col == NIL(sm_element)) {
  287.         sm_delrow(A, prow->row_num);
  288.         }
  289.     }
  290.  
  291.     /* discard the column -- we already threw away the elements */ 
  292.     A->cols[i] = NIL(sm_col);
  293.     dll_unlink(pcol, A->first_col, A->last_col, 
  294.                 next_col, prev_col, A->ncols);
  295.     pcol->first_row = pcol->last_row = NIL(sm_element);
  296.     sm_col_free(pcol);
  297.     }
  298. }
  299.  
  300. void sm_copy_row(register sm_matrix *dest, int dest_row, sm_row *prow)
  301. {
  302.     register sm_element *p;
  303.  
  304.     for(p = prow->first_col; p != 0; p = p->next_col) {
  305.     (void) sm_insert(dest, dest_row, p->col_num);
  306.     }
  307. }
  308.  
  309.  
  310. void sm_copy_col(register sm_matrix *dest, int dest_col, sm_col *pcol)
  311. {
  312.     register sm_element *p;
  313.  
  314.     for(p = pcol->first_row; p != 0; p = p->next_row) {
  315.     (void) sm_insert(dest, dest_col, p->row_num);
  316.     }
  317. }
  318.  
  319.  
  320. sm_row *sm_longest_row(sm_matrix *A)
  321. {
  322.     register sm_row *large_row, *prow;
  323.     register int max_length;
  324.  
  325.     max_length = 0;
  326.     large_row = NIL(sm_row);
  327.     for(prow = A->first_row; prow != 0; prow = prow->next_row) {
  328.     if (prow->length > max_length) {
  329.         max_length = prow->length;
  330.         large_row = prow;
  331.     }
  332.     }
  333.     return large_row;
  334. }
  335.  
  336.  
  337. sm_col *sm_longest_col(sm_matrix *A)
  338. {
  339.     register sm_col *large_col, *pcol;
  340.     register int max_length;
  341.  
  342.     max_length = 0;
  343.     large_col = NIL(sm_col);
  344.     for(pcol = A->first_col; pcol != 0; pcol = pcol->next_col) {
  345.     if (pcol->length > max_length) {
  346.         max_length = pcol->length;
  347.         large_col = pcol;
  348.     }
  349.     }
  350.     return large_col;
  351. }
  352.  
  353.  
  354. int sm_num_elements(sm_matrix *A)
  355. {
  356.     register sm_row *prow;
  357.     register int num;
  358.  
  359.     num = 0;
  360.     sm_foreach_row(A, prow) {
  361.     num += prow->length;
  362.     }
  363.     return num;
  364. }
  365.  
  366. int sm_read(FILE *fp, sm_matrix **A)
  367. {
  368.     int i, j, err;
  369.  
  370.     *A = sm_alloc();
  371.     while (! feof(fp)) {
  372.     err = fscanf(fp, "%d %d", &i, &j);
  373.     if (err == EOF) {
  374.         return 1;
  375.     } else if (err != 2) {
  376.         return 0;
  377.     }
  378.     (void) sm_insert(*A, i, j);
  379.     }
  380.     return 1;
  381. }
  382.  
  383.  
  384. int 
  385. sm_read_compressed(FILE *fp, sm_matrix **A)
  386. {
  387.     int i, j, k, nrows, ncols;
  388.     unsigned long x;
  389.  
  390.     *A = sm_alloc();
  391.     if (fscanf(fp, "%d %d", &nrows, &ncols) != 2) {
  392.     return 0;
  393.     }
  394.     sm_resize(*A, nrows, ncols);
  395.  
  396.     for(i = 0; i < nrows; i++) {
  397.     if (fscanf(fp, "%lx", &x) != 1) {
  398.         return 0;
  399.     }
  400.     for(j = 0; j < ncols; j += 32) {
  401.         if (fscanf(fp, "%lx", &x) != 1) {
  402.         return 0;
  403.         }
  404.         for(k = j; x != 0; x >>= 1, k++) {
  405.         if (x & 1) {
  406.             (void) sm_insert(*A, i, k);
  407.         }
  408.         }
  409.     }
  410.     }
  411.     return 1;
  412. }
  413.  
  414.  
  415. void sm_write(FILE *fp, sm_matrix *A)
  416. {
  417.     register sm_row *prow;
  418.     register sm_element *p;
  419.  
  420.     for(prow = A->first_row; prow != 0; prow = prow->next_row) {
  421.     for(p = prow->first_col; p != 0; p = p->next_col) {
  422.         (void) fprintf(fp, "%d %d\n", p->row_num, p->col_num);
  423.     }
  424.     }
  425. }
  426.  
  427. void sm_print(FILE *fp, sm_matrix *A)
  428. {
  429.     register sm_row *prow;
  430.     register sm_col *pcol;
  431.     int c;
  432.  
  433.     if (A->last_col->col_num >= 100) {
  434.     (void) fprintf(fp, "    ");
  435.     for(pcol = A->first_col; pcol != 0; pcol = pcol->next_col) {
  436.         (void) fprintf(fp, "%d", (pcol->col_num / 100)%10);
  437.     }
  438.     putc('\n', fp);
  439.     }
  440.  
  441.     if (A->last_col->col_num >= 10) {
  442.     (void) fprintf(fp, "    ");
  443.     for(pcol = A->first_col; pcol != 0; pcol = pcol->next_col) {
  444.         (void) fprintf(fp, "%d", (pcol->col_num / 10)%10);
  445.     }
  446.     putc('\n', fp);
  447.     }
  448.  
  449.     (void) fprintf(fp, "    ");
  450.     for(pcol = A->first_col; pcol != 0; pcol = pcol->next_col) {
  451.     (void) fprintf(fp, "%d", pcol->col_num % 10);
  452.     }
  453.     putc('\n', fp);
  454.  
  455.     (void) fprintf(fp, "    ");
  456.     for(pcol = A->first_col; pcol != 0; pcol = pcol->next_col) {
  457.     (void) fprintf(fp, "-");
  458.     }
  459.     putc('\n', fp);
  460.  
  461.     for(prow = A->first_row; prow != 0; prow = prow->next_row) {
  462.     (void) fprintf(fp, "%3d:", prow->row_num);
  463.  
  464.     for(pcol = A->first_col; pcol != 0; pcol = pcol->next_col) {
  465.         c = sm_row_find(prow, pcol->col_num) ? '1' : '.';
  466.         putc(c, fp);
  467.     }
  468.     putc('\n', fp);
  469.     }
  470. }
  471.  
  472.  
  473. void sm_dump(sm_matrix *A, char *s, int max)
  474. {
  475.     FILE *fp = stdout;
  476.  
  477.     (void) fprintf(fp, "%s %d rows by %d cols\n", s, A->nrows, A->ncols);
  478.     if (A->nrows < max) {
  479.     sm_print(fp, A);
  480.     }
  481. }
  482.  
  483. void
  484. sm_cleanup()
  485. {
  486. #ifdef FAST_AND_LOOSE
  487.     register sm_element *p, *pnext;
  488.     register sm_row *prow, *pnextrow;
  489.     register sm_col *pcol, *pnextcol;
  490.  
  491.     for(p = sm_element_freelist; p != 0; p = pnext) {
  492.     pnext = p->next_col;
  493.     FREE(p);
  494.     }
  495.     sm_element_freelist = 0;
  496.  
  497.     for(prow = sm_row_freelist; prow != 0; prow = pnextrow) {
  498.     pnextrow = prow->next_row;
  499.     FREE(prow);
  500.     }
  501.     sm_row_freelist = 0;
  502.  
  503.     for(pcol = sm_col_freelist; pcol != 0; pcol = pnextcol) {
  504.     pnextcol = pcol->next_col;
  505.     FREE(pcol);
  506.     }
  507.     sm_col_freelist = 0;
  508. #endif
  509. }
  510.